#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <glpk.h>
using namespace std;

vector<double> Find_XI(vector<double> _beta, int i, double eps)
{
    vector<double> X;
    for(double j = _beta[i];j < _beta[i+1];j+=eps)
    {
        X.push_back(j);
    }
    return X;
}
/*
double fraction_rho(vector<double> _r, int _Up, int _Down)
{
    double sum = 1.0;
    for(int j=_Down;j<=_Up;j++)
    {
        double tmp1 = 0.0;
        for(int z=1;z<=_Up;z++)
        {
            tmp1 += _r[z];
        }
        double tmp2 = 0.0;
        for(int z=1;z<_Down;z++)
        {
            tmp2 += _r[z];
        }
        sum *= (_r[j] / (tmp1 - tmp2));
        return sum;
    }
}

double rho_leq(vector<double> _r, int index)
{
    double sum = 0.0;
    for(int i=1;i<=index;i++)
    {
        sum += _r[i];
    }
    return sum;
}
*/
int main()
{
	double final_ans = 0.0;
	int index = 0, final_index = index;

	for(double rho1=0.6831;rho1<=0.6831;rho1+=0.001)
	{
		for(double rho2=0.2342;rho2<=0.2342;rho2+=0.001)
		{
			for(double rho3=0.0711;rho3<=0.0711 && ((rho1 + rho2 + rho3) <= 1.0);rho3+=0.001)
			{
				for(double c=0.6906;c<=0.6906;c+=0.001)
				{
					// double rho4 = 1.0 - rho1 - rho2 - rho3;
					ofstream file;
					file.open("LP.txt", ios::trunc);

					/**Param**/
					file << "param k;" << endl;
					file << "param rho1;" << endl;
					file << "param rho2;" << endl;
					file << "param rho3;" << endl;
					file << "param rho4;" << endl;
					file << "set I:= {0 .. k};" << endl;
					file << "set II:= {1 .. k};" << endl;
					file << "set III:= {0 .. k-1};" << endl;
					file << "param beta{I};" << endl;
					file << "param exp_beta_rho1{I};" << endl; // exp(-beta[i] * rho1)
					file << "param exp_beta_rho2{I};" << endl; // exp(-beta[i] * rho2)
					file << "param exp_beta_rho3{I};" << endl; // exp(-beta[i] * rho3)


					file << "param exp_beta_rho4{I};" << endl; // exp(-beta[i] * rho4)
					file << "param exp_beta_rho{I};" << endl; // exp(-beta[i] * (rho1+rho2))
					file << "param exp_beta_rho123{I};" << endl; // exp(-beta[i] * (rho1 + rho2 + rho3))
					file << "param exp_beta_rho23{I};" << endl; // exp(-beta[i] * (rho2 + rho3))
					int k = 2500;
					double rho4 = 1.0 - rho1 - rho2 - rho3;
					// double rho1 = 0.6828, rho2 = 0.2335, rho3 = 0.0695 , rho4 = 1.0 - rho1 - rho2 - rho3, c = 0.6910;
					vector<double> beta(k+1);
					for(int i=0;i<=k;i++)
					{
						beta[i] = c + 0.01 * i;
					}


					/**Var**/
					file << "var resultoffourphase;" << endl;
					file << "var Lambda1;" << endl;
					file << "var Lambda2;" << endl;
					file << "var Lambda3;" << endl;
					file << "var Lambda4;" << endl;
					file << "var Delta{I};" << endl;
					file << "var delta{III};" << endl;
					file << "var v{I};" << endl;

					/**Objective Function**/
					file << "minimize z: resultoffourphase;" << endl;
					
					/**Constraint**/
					file << "subject to result4: resultoffourphase >= Lambda1 + Lambda2 + Lambda3 + Lambda4;" << endl;
					file << "subject to Low_Lambda1: " << 1.0 - exp(-beta[0] * rho1) << " * (v[0] + Delta[0] / beta[0]) <= Lambda1;" << endl;
					file << "subject to Low_Lambda2: " << "sum {i in II} (exp_beta_rho1[i-1] - exp_beta_rho1[i]) * (1 - exp_beta_rho2[i]) * (beta[i-1] / beta[i] * v[i] + Delta[i-1]/beta[i]) <= Lambda2;" << endl;
					file << "subject to Low_Lambda3: " << "sum {i in II} (1 - exp_beta_rho3[i]) * (beta[i-1] / beta[i] * v[i] + Delta[i-1]/beta[i]) * (" << exp(-rho1 * c) << "*(exp_beta_rho2[i-1] - exp_beta_rho2[i]) - rho2 / (rho1 + rho2) * (exp_beta_rho[i-1] - exp_beta_rho[i])) <= Lambda3;" << endl; 
					file << "subject to Low_Lambda4: " << "sum {i in II} (((rho1 / (rho1 + rho2)) * (exp_beta_rho[0] * (exp_beta_rho3[i-1] - exp_beta_rho3[i]) - (exp_beta_rho123[i-1] - exp_beta_rho123[i])) + (rho2 / (rho2 + rho3)) * (exp_beta_rho1[0] * (exp_beta_rho23[i-1] - exp_beta_rho23[i]) - (exp_beta_rho123[i-1] - exp_beta_rho123[i])) + (rho1 / (rho1 + rho2 + rho3)) * (rho2 / (rho2 + rho3)) * (exp_beta_rho123[i-1] - exp_beta_rho123[i])) * (1 - exp_beta_rho4[i]) * (beta[i-1] / beta[i] * v[i] + Delta[i-1]/beta[i]))" << " <= Lambda4;" << endl;
					
					file << "subject to OPT: " << "1 <= v[k] + Delta[0] + sum {i in III} delta[i];" << endl;
					for(int i=0;i<k;i++)
					{
						file << "subject to Up_delta_one_C" << i << ": delta[" << i << "] <= " << 1.0 - exp(-beta[i+1]) << " * (v[" << i << "] - v[" << i+1 << "]);" << endl;
					}
					for(int i=0;i<k;i++)
					{
						file << "subject to Up_delta_two_C" << i << ": delta[" << i << "] <= Delta[" << i+1 << "] - Delta[" << i << "];" << endl; 
					}
					for(int i=0;i<k;i++)
					{
						vector<double> Xi = Find_XI(beta, i, 0.001);
						int index1 = 0;
						for(auto it=Xi.begin();it!=Xi.end();it++)
						{
							double gamma = (*it - beta[i]) / (exp(-*it) - exp(-beta[i+1]));
							file << "subject to Up_delta_three_C" << i << "_" << index1 << ": " << gamma << " * delta[" << i << "] <= Delta[" << i+1 << "] - Delta[" << i << "] - " << beta[i] - gamma * (1.0 - exp(-*it)) << " * (v[" << i << "] - v[" << i+1 << "]);" << endl;
							index1++;
						}
					}
					for(int i=0;i<k;i++)
					{
						file << "subject to V_C" << i << ": v[" << i+1 << "] <= v[" << i << "];" << endl;
					}
					file << "subject to V_C" << k << ": v[" << k << "] >= 0;" << endl;
					file << "subject to Low_Delta_C0: Delta[0] >= 0;" << endl;
					for(int i=0;i<k;i++)
					{
						file << "subject to Low_delta_C" << i << ": delta[" << i << "] >= 0;" << endl;
					}
					/***Data*/
					file << "data;" << endl;
					file << "param k := " << k << ";" << endl;
					file << "param rho1 := " << rho1 << ";" << endl;
					file << "param rho2 := " << rho2 << ";" << endl;
					file << "param rho3 := " << rho3 << ";" << endl;
					file << "param rho4 := " << rho4 << ";" << endl;

					file << "param beta :=" << endl;
					for(int i=0;i<=k;i++)
					{
						file << i << " " << beta[i] << endl;
					}
					file << ";" << endl;

					file << "param exp_beta_rho1 :=" << endl;
					for(int i=0;i<=k;i++)
					{
						file << i << " " << exp(-beta[i] * rho1) << endl;
					}
					file << ";" << endl;

					file << "param exp_beta_rho2 :=" << endl;
					for(int i=0;i<=k;i++)
					{
						file << i << " " << exp(-beta[i] * rho2) << endl;
					}
					file << ";" << endl;

					file << "param exp_beta_rho3 :=" << endl;
					for(int i=0;i<=k;i++)
					{
						file << i << " " << exp(-beta[i] * rho3) << endl;
					}
					file << ";" << endl;

					file << "param exp_beta_rho4 :=" << endl;
					for(int i=0;i<=k;i++)
					{
						file << i << " " << exp(-beta[i] * rho4) << endl;
					}
					file << ";" << endl;

					// for(int t=3;t<=m;t++)
					// {
					//     file << "param exp_beta_rho" << t << " :=" << endl;
					//     for(int i=0;i<=k;i++)
					//     {
					//         file << i << " " << exp(-beta[i] * rho[t]) << endl;
					//     }
					//     file << ";" << endl;
					// }

					file << "param exp_beta_rho :=" << endl;
					for(int i=0;i<=k;i++)
					{
						file << i << " " << exp(-beta[i] * (rho1+rho2)) << endl;
					}
					file << ";" << endl;

					file << "param exp_beta_rho123 :=" << endl;
					for(int i=0;i<=k;i++)
					{
						file << i << " " << exp(-beta[i] * (rho1+rho2+rho3)) << endl;
					}
					file << ";" << endl;

					file << "param exp_beta_rho23 :=" << endl;
					for(int i=0;i<=k;i++)
					{
						file << i << " " << exp(-beta[i] * (rho2 + rho3)) << endl;
					}
					file << ";" << endl;
					file.close();
					glp_prob *mip;
					glp_tran *tran;
					int ret;
					mip = glp_create_prob();
					tran = glp_mpl_alloc_wksp();
					ret = glp_mpl_read_model(tran, "LP.txt", 0);
					if(ret != 0)
					{
						fprintf(stderr, "Error on translating model\n");
						// goto skip;
					}
					ret = glp_mpl_generate(tran, NULL);
					if(ret != 0)
					{
						fprintf(stderr, "Error on generating model\n");
						// goto skip;
					}
					glp_mpl_build_prob(tran, mip);
					glp_scale_prob(mip, GLP_SF_AUTO);
					glp_smcp parm;
					glp_init_smcp(&parm);
					parm.meth = GLP_DUALP;
					// parm.meth = GLP_PRIMAL;
					parm.presolve = GLP_ON;
					glp_simplex(mip, &parm);
					// glp_simplex(mip, NULL);
					// glp_intopt(mip, NULL);
					ret = glp_mpl_postsolve(tran, mip, GLP_SOL);
					if(ret != 0)
					{
						fprintf(stderr, "Error on postsolving model\n");
					}
					int TMPTMP = glp_print_sol(mip, "123Result.txt");
					// ret = glp_print_sol(mip, NULL);
					double ans = glp_get_obj_val(mip);
					final_ans = max(final_ans, ans);
					file.open("Result.txt", ios::app);
					file << "rho1: " << rho1 << ", rho2: " << rho2 << ", rho3: " << rho3 << ", c: " << c << ", ratio: " << ans << endl;
					file.close();
					index++;
					if(final_ans == ans)
					{
						final_index = index;
					}
					glp_mpl_free_wksp(tran);
					glp_delete_prob(mip);
			/**	
				skip: glp_mpl_free_wksp(tran);
					  glp_delete_prob(mip);	
			**/

	
				}
			}
		}
	}
	cout << "The Maximum ratio is " << final_ans << ", and the index is " << final_index << endl;
	return 0;
}
